Carta EWMA

Parâmetros do modelo

nH0 <- 100
nH1 <- 200
phi_parametro <- 0.2

coeficientes <- function(phi) {
  # βARMA: the model from Rocha and Cribari-Neto (2009, 2017)
  #        is obtained by setting coefs$d = 0 and d = FALSE and error.scale = 1 (predictive scale)
  list(alpha = 0, beta = NULL, phi = phi, theta = NULL, d = 0, nu = 20)
}

barma.sim <- function(n, phi, seed, y.start = NULL){
  BARFIMA.sim(
    n = n,
    coefs = coeficientes(phi),
    y.start = y.start,
    error.scale = 1,
    complete = F,
    seed = seed
  )
}

barma.phi_estimado <- function(yt, alpha = 0, nu = 20, phi = 0.1){
  BARFIMA.fit(
    yt = yt,
    start = list(alpha = alpha, nu = nu, phi = phi),
    p = 1,
    d = FALSE,
    error.scale = 1,
    report = F
  )$coefficients["phi"][[1]]
}

barma.residuos <- function(yt, phi_estimado){
  BARFIMA.extract(
    yt = yt,
    coefs = coeficientes(phi_estimado),
    llk = F,
    info = F,
    error.scale = 1
  )$error
}
# Função para forçar a execução de `BARFIMA.extract` até que ela funcione
so_quero_que_funcione <- function(func, debug = T, max = 1000, ...) {
  FINALMENTE <- F
  contador <- 0
  
  while (!FINALMENTE) {
    contador <- contador + 1
    if (debug) print(paste("Tentativa:", contador))
    
    resultado <- tryCatch({
      func(...)
    }, error = function(err) {
      if (contador >= max) {
        stop("Deu ruim, número máximo de tentativas atingido")
      }
      
      if (any(grepl("\\.btsr\\.extract\\(", deparse(err$trace$call)))) {
        # Ignora erro de extração de resíduos
        return(NULL)
      }
      
      stop(err)
    })
    
    if (!is.null(resultado)) {
      FINALMENTE <- T
    }
  }
  
  return(resultado)
}
ewma_qcc <- function(amostra_inicial, lambda, amostra_subsequente) {
  ew <- ewma(amostra_inicial, lambda = lambda, nsigmas = 3, plot = F, newdata = amostra_subsequente)
  
  registros <- (nH0 + 1):(nH0 + nH1)
  list(
    ewma = as.numeric(ew$y[registros]),
    LI = ew$limits[, 1][registros],
    LS = ew$limits[, 2][registros]
  )
}

Simulação de amostras subsequentes

gera_amostras_subsequentes <- function() {
  lambda <- 0.2
  
  amostra_controle <- barma.sim(
    n = nH0,
    phi = phi_parametro,
    seed = 1337
  )
  
  ultima_observacao <- amostra_controle[nH0]
  phi_estimado <- barma.phi_estimado(amostra_controle)
  residuos_controle <- barma.residuos(amostra_controle, phi_estimado)
  
  data.frame(
      novo_phi = seq(0, 0.8, by = 0.1)
    ) %>%
    rowwise() %>%
    mutate(
      observacao = list(1:nH1),
      dados = list(barma.sim(
        n = nH1,
        phi = novo_phi,
        seed = 1337,
        y.start = ultima_observacao
      )),
      residuos = list(
        barma.residuos(dados, phi_estimado)
      ),
      qcc = list(ewma_qcc(residuos_controle, lambda, residuos)),
      ewma = list(qcc$ewma),
      LI = list(qcc$LI),
      LS = list(qcc$LS),
      fora_de_controle = list(ewma < LI | ewma > LS),
      total_fora_de_controle = sum(fora_de_controle),
      fracao_fora_de_controle = total_fora_de_controle / nH1
    )
}


amostras_subsequentes <- so_quero_que_funcione(gera_amostras_subsequentes)
## [1] "Tentativa: 1"
amostras_subsequentes %>%
  group_by(novo_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    .groups = "drop"
  )

Cartas

ggplotly(
  amostras_subsequentes %>%
    unnest(
      cols = c(novo_phi, observacao, ewma, LI, LS)
    ) %>%
    ggplot() +
    geom_line(aes(x = observacao, y = LI, color = "Limite Inferior"), linetype = "dashed") +
    geom_line(aes(x = observacao, y = LS, color = "Limite Superior"), linetype = "dashed") +
    geom_line(aes(x = observacao, y = ewma, color = "EWMA")) +
    geom_point(data = . %>% filter(ewma < LI | ewma > LS),
               aes(x = observacao, y = ewma, color = "Fora de controle"), size = 1.5, shape = 4) +
    labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
    scale_color_manual(
      values = c(
        "Limite Inferior" = adjustcolor("red", alpha.f = 0.5),
        "Limite Superior" = adjustcolor("red", alpha.f = 0.5),
        "EWMA" = adjustcolor("black", alpha.f = 0.8),
        "Fora de controle" = adjustcolor("red", alpha.f = 0.4)
      )
    ) +
    labs(
      title = "EWMA para resíduos βAR(1), com λ = 0.2",
    ) +
    theme.base +
    theme(
      legend.position = "bottom",
      legend.title = element_blank()
    ) +
    facet_wrap(~novo_phi)
) %>%
  plotly.base

Simulação de Monte Carlo

amostras_monte_carlo_gerador <- function(){
  numero_de_execucoes <- 200
  
  expand.grid(
    k = 1:numero_de_execucoes,
    novo_phi = seq(0.2, 0.6, by = 0.1),
    lambda = seq(0.1, 0.4, by = 0.1)
  ) %>%
  mutate(
    id = row_number()
  ) %>%
  rowwise() %>%
  mutate(
    amostra_controle = list(barma.sim(
      n = nH0,
      phi = phi_parametro,
      seed = id
    )),
    ultima_observacao = amostra_controle[nH0],
    phi_estimado = barma.phi_estimado(amostra_controle),
    residuos_controle = list(
      barma.residuos(amostra_controle, phi_estimado)
    ),
    dados = list(barma.sim(
      n = nH1,
      phi = novo_phi,
      y.start = ultima_observacao,
      seed = id + 1337E4
    )),
    residuos = list(
      barma.residuos(dados, phi_estimado)
    ),
    qcc = list(ewma_qcc(residuos_controle, lambda, residuos)),
    ewma = list(qcc$ewma),
    LI = list(qcc$LI),
    LS = list(qcc$LS),
    fora_de_controle = list(ewma < LI | ewma > LS),
    total_fora_de_controle = sum(fora_de_controle, na.rm = TRUE),
    fracao_fora_de_controle = total_fora_de_controle / nH1
  )
}

if (file.exists("cache_simulacao_ewma.RData")) {
  amostras_monte_carlo <- readRDS("cache_simulacao_ewma.RData")
} else {
  amostras_monte_carlo <- so_quero_que_funcione(amostras_monte_carlo_gerador)
  saveRDS(amostras_monte_carlo, file = "cache_simulacao_ewma.RData")
}

Resumo

amostras_monte_carlo_resumo <- amostras_monte_carlo %>%
  group_by(lambda, novo_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    min = min(fracao_fora_de_controle),
    max = max(fracao_fora_de_controle),
    .groups = "drop"
  )

amostras_monte_carlo_resumo

Gráficos

ggplotly(
  amostras_monte_carlo_resumo %>%
      ggplot(aes(x = novo_phi, y = mean, color = factor(lambda))) +
      geom_line() +
      geom_point() +
      labs(
        x = "Valores de Φ",
        y = "Fração de pontos fora de controle",
        color = "Valores de λ",
        title = "Fração de pontos fora de controle por Φ e λ"
      ) +
      theme.base
) %>%
  plotly.base

Box plot

ggplotly(
  amostras_monte_carlo %>%
    ggplot(aes(x = factor(novo_phi), y = fracao_fora_de_controle, fill = factor(lambda))) +
    geom_boxplot() +
    labs(
      x = "Valores de Φ",
      y = "Fração de pontos fora de controle",
      fill = "Valores de λ",
        title = "Fração de pontos fora de controle por Φ e λ"
    ) +
    theme.base
  ) %>%
  plotly.base %>%
  layout(boxmode = 'group')